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ABSTRACT 

In this paper we study the feathering substructures along spiral arms by considering the perturba- 
tional gas response to a spiral shock. Feathers are density fluctuations that jut out from the spiral 
arm to the inter-arm region at pitch angles given by the quantum numbers of the doubly-periodic 
structure. In a localized asymptotic approximation, related to the shearing sheet except that the 
inhomogeneities occur in space rather than in time, we derive the linearized perturbation equations 
for a razor-thin disk with turbulent interstellar gas, frozen-in magnetic field, and gaseous self-gravity. 
Apart from the modal quantum numbers, the individual normal modes of the system depend on seven 
dimensionless quantities that characterize the underlying time-independent axisymmetric state plus 
its steady, nonlinear, two-armed spiral-shock (TASS) response to a hypothesized background density- 
wave supported by the disk stars of the galaxy. We show that some of these normal modes have 
positive growth rates. Their over-density contours in the post-shock region are very reminiscent of 
observed feathering substructures in full magnetohydrodynamic (MHD) simulations. The feathering 
substructures are parasitic instabilities intrinsic to the system; thus, their study not only provides 
potential diagnostics for important parameters that characterize the interstellar medium of external 
galaxies, but also yields a deeper understanding of the basic mechanism that drives the formation of 
the giant molecular clouds (GMCs) and the OB stars that outline observed grand-design spirals. 
Subject headings: Galaxies: ISM, Galaxies: Structure, Instabilities, ISM: Kinematics and Dynamics, 
ISM: Magnetic Fields, Magnetohydrodynamics: MHD 



1. INTRODUCTION 

Spiral structures in nearby galaxies have fascinated 
astronomers since Lord Rosse's observations of M51 in 
1845. The underpinning for a theoretical understand- 
ing of the phenomenon in terms of density waves has 
existed about 50 years (see iLin fe SLvul 119641 and refer- 
ences therein). Improved imaging technology and tech- 
niques reveal many substructures associated with the spi- 
ral arms. Here, we focus on the quasi-regularly spaced 
density fluctu ations ide ntified in the litera ture as feathers 
(ILvndsl [T970 N ) or spurs (Elmcgrccn 1980). Observation- 
ally, the feathers are extinction substructure commonly 
found in the optical band im ages among spiral galax- 
ies, e.g., lLa Vigne et "all (l200l. For example, a Hubble 
Heritage image of M51 (IScoville fc Rectorl l20f)l . shows 
many feathers (i.e., darkened dust lane in the optical im- 
age) projected into the inter-arm region from the primary 
dust lane. There are also examples showing the feathers 
in infrared (e.g., the 8/im image of M81 from Spitzers 
Space Telescope) or sub-millimeter wavele ngths, such as 
detec tion of CO emission in M51 feathers (jCorder et al.1 
2008). Therefore, the relationship between the feathers 
and the underlying interstellar medium (molecular and 
atomic gas, dust, magnetic field, etc.) may hold the key 
to an understanding of the formation of the GMCs and 
OB stars that delineate the arms of spiral galaxies. 

There are two points of view regarding the background 
structures of spiral galaxies. The first is the hypothesis of 
quasi-stationary spiral structure (QSSS) that attributes 
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the origin of spiral structure to the normal modes of 
the disk stars of a flattened galaxy. This point of view 
seems consistent with the observational finding that spi- 
ral galaxies, which look fragmentary, multi-armed, and 
even flocculent at optical or blue wavelengths, never- 
theless have, in 2.1/im images, the grand-design two- 
armed spiral stru cture (TASS) that underl i es the QSSS 
hypo t hesis (e.g., iBlock fc Wainscoatl 119911: IBlock et al.l 
119941: IBlock. Elmegreen. fc Wainscoatl 119961 ) . The sec- 
ond comes from numerical simulations that show non- 
linear effect s saturating the growth of unstable normal 
modes (e.g., Scllwood 2012) lead to spiral patterns that 
are locally transient. These are dichotomies of long 
standing that we do not address in the present paper, 
which focuses on the substructures that arise from the re- 
sponse of the self-gravitating and magnetized interstellar 
medium even to the steady forcing associated with the 
classic QSSS hypothesis. 

Theoretical understanding of the substructures is 
also confused. Explanations encompass both irregu- 
lar causes such as swing- amplified shea ring instabili- 
ties (e.g., IGoldreich fc Lvnden-Belll Il965l ). and regular 
causes such as gravitation al instabilities (e.g., [Balbus 
ll988tlKim fc Ostrikerll2002D initiated b y a TASS pattern 
(jRobertsI 119691: iRoberts fc Yuanlll970D . Another possi- 
bility is that spurs arise as response of the disk stars 
to ov e r-dense regions like G MCs (e.g.. IJulian fc Toomrd 
119661: ID'Onghia et al.l I2012D . The last possibility will 
lead, however, to spurs with characteristic inclinations 
that co-rotate with the local material velocity of the 
GMCS, which would not have an obvious correlation with 



2 



Lee and Shu 



the spiral pattern of the older disk stars. Elmegreen's 
(1980) conclusion that spurs have characteristic incli- 
nations that correlate with / band images of the older 
disk stellar population suggests that feathering is best 
described as a long-lived phenomenon, intimately con- 
nected with the underlying spiral structure of disk galax- 
ies. We adopt this hypothesis for the analysis of the 
present paper, and do not speculate on the changes nec- 
essary if spiral patterns are short-lived with spiral arms 
persistent only in a statistical sense. 

The QSSS hypothesis implicitly underlies many nu- 
merical si mulations in rec e nt y e ars on this sub- 
ject (e.g.. iKim fc Ostrikerl l2002t iChakrabarti et all 
2001 Kim fc Ostrikerl 12001 [Shetty fc Ostrikerl 120061: 
Dobbs fc Bonnelll 120061 |Dobbs] |2008|) . These sophisti- 
cated simulations include MHD, self-gravity, ISM phases, 
etc. They provide a detailed time evolution of how GMCs 
can be formed by the fragmentation and agglomeration 
of interstellar gas by local Jeans instability. However, 
due to computational limitations, the behavior of the 
system is followed for only a few orbital times. Also, 
as we shall see, given that seven dimensionless numbers 
form the irreducible set that characterizes the instabil- 
ity of the system, a comprehensive survey of parameter 
space by numerical simulation is clearly out of the ques- 
tion for the foreseeable future. On the other hand, most 
theoretical linearized-stability analyses along the same 
line of thought include restrictive assumptions such as 
an arbitrary background profile, and/or a shearing-sheet 
approach, and/or a lack of gaseous self-gravity and/or 
magnetic fields (e.g... iDwarkadas fc Baibuslll996l ). These 
simplifications compromise the applicability of the analy- 
sis if we wish ultimately to use the theory as a diagnostic 
of the physical conditions in real systems. Our aim here 
is to rectify these shortcomings. 

In this paper we formulate and solve the basic equa- 
tions that govern the formation of feathers through the 
instability of a galactic spiral shock when the roles 
of gaseous self-gravity and magnetic fi eld are i n clude d 
within the original TASS framework of IRobertsI (|1969h . 
We work in the frame that corotates with pattern speed 
f2 p of the spiral gravitational field of the background stel- 
lar disk, in which the TASS pattern is independent of 
time t and asymptotically one-dimensional (i.e., varia- 
tions only in the direction perpendicular to spiral arms) . 
By transforming the governing nonlinear equations to a 
spiral coordinate system (77, £) with 77 varying perpendic- 
ular to spiral arms, and £ along them, we write down the 
asymptotic equations that govern nonlinear behavior in 
which the underlying TASS pattern varies only in 77 but 
the parasitic perturbations above the TASS pattern can 
vary in all three variables (77, £, t). The self-consistency of 
the asymptotic approximation then requires us to impose 
that single-valued perturbations are doubly periodic in 
(77, £) when we linearize in the amplitude of the perturba- 
tions relative to the TASS state. This double-periodicity 
is characterized by two integers (quantum numbers): m 
= the number of stellar spiral arms in a complete cir- 
cle around the galaxy with m assumed to equal to 2 in 
practice, and / = the number of feathers as we go along 
a spiral arm that would take us to the next spiral arm 
(half-way circumferentially around the galaxy if m = 2) 
if we were to go instead in the direction perpendicular to 
a spiral arm. 



Our calculation on the TASS part of the problem dif- 
fers from the original Roberts work in that we include 
frozen- in magnetic fields (as did lRoberts fc Yuanlll970l) 
and gaseous self-gravity (Ostriker and Kim's 2002 anal- 
ysis included only the self-gravity of the feathering per- 
turbations, and not its effect on the underlying TASS 
state) . When we also include the effect of turbulent mo- 
tions of the interstellar gas, modelled as a "logatropic" 
gas (pressure P proportional to the logarithm of the den- 
sity p), there are seven dimensionless, irreducible, num- 
bers that characterize the TASS state: (1) the ratio of 
the circular frequency to the epicyclic frequency O/k; 
(2) the sine of the inclination of the stellar spiral arms, 
sini; (3) the dimensionless Doppler-shifted frequency at 
which gas rotating at its circular angular speed f2 meets 
m stellar spiral arms that each rotate at angular speed 
ftp, —v = 771 (f2 — £I p )/k; (4) the amplitude of the stellar 
spiral gravitational field as a fraction of the axisymmetric 
radial gravitational field, F; (5) the dimensionless mea- 
sure of the gas surface density, a; (6) the dimensionless 
measure of the mean gas turbulent speed, Xt] and (7) the 
dimensionless measure of the mean Alfven speed of the 
magnetized interstellar medium, xa- 

The plan of the paper is as follows. In <|2j we write 
down the basic MHD equations in the spiral coordinates. 
In we obtain by the shooting method the 1-D nonlin- 
ear TASS solution in 77 (across the spiral arm) modified 
from the Roberts-style analysis by the inclusion of mag- 
netic fields and gaseous self-gravity. In SjH we derive 
the equations that govern linearized, time-dependent, 2- 
D perturbations on top of the background TASS pattern. 
Because the basic reference state depends only on 77 and 
not £ nor t, the linearized perturbations can be taken 
to be oscillatory (with complex frequency ur + iLO\) in 
time t and (with real dimensionless wavenumber I) in the 
spatial dimension £, but with dependences on the spatial 
dimension 77 that satisfy ordinary differential equations. 
Generality requires us to consider oscillatory perturba- 
tions in the position of the shock. When the appropriate 
jump conditions (due to the corrugation of the shock 
front) are imposed on top of the condition of double- 
periodicity, we obtain the real and imaginary parts of 
the perturbation frequency, wr and u>\, as eigenvalues of 
the problem when the quantum numbers m and / are 
specified, together with the numerical values of Q/k, is, 
sini, F, a, Xt, and xa- In $5] we give a sample result. In 
£j6] we discuss the physical meaning of the result and give 
our conclusions. 

2. BASIC EQUATIONS AND GEOMETRY 

We first write down the basic equations for the problem 
from the two-dimensional, time-dependent, ideal MHD 
equations in the rotating frame. We identify the ax- 
isymmetric, time-independent solution as the zeroth or- 
der state. We then introduce the tight-winding spiral 
arm approximation and obtain the I s * order (in term 
of sin i) nonlinear TASS state in the sense of IRobertsI 
(fl969l) . modified for gaseous self-gravity and the pres- 
ence of frozen-in magnetic fields. 

2.1. Basic Equations 

In cylindrical polar coordinates (tn, ip, z), we denote S, 
7i ro and u v as, respectively, the gas surface density, w- 
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and (^-components of the fluid velocity in a razor-thin 
flat disk. The continuity and momentum equations, in 
a rotating frame with angular rate of the spiral pattern, 



J7p, can be written as, 



dt 



Id Id 

- — H — (£%,) = 0; 

vj dvj vj dip 



dUr- 



+ 



dUr; 



dvj vj dip 
dV cS 



v 
w 



_ldU 

™ E dvo dvo 



(i) 



(2) 



fraction of the axisymmetric gravitational acceleration: 



F 



\k„\A 
vjQ? ' 



(8) 



A typical number quoted in the literature is F — 
5 to 10%. Although infrared images may indicate 
stronger fractions compared to the background stel- 
lar disk, especially in the outer disk, it should be re- 
called that the denominator in the definition of F in- 
cludes the force contribu t ion fro m the dark matter halo. 
iShu. Milione. fc Robertsl (| 19731 ) (hereafter SMR) show, 
however, that the real measure of nonlinearity of the 
gaseous forcing is given by the combination, 



du v 
~~dt 



du v 
dvo 

i an 

vjYj dp 



du v u^u 



VJ U'tf 



vj d(p 

i aveff 

vj dp 



(3) 



where II is the vertically- integrated gas pressure, Fru,Fip 
are the two horizontal components of the Lorentz force 
per unit mass. The last terms in equations ([2]) and ([3]) 
are the Coriolis accelerations associated with being in a 
frame of reference that rotates at angular speed f2 p . We 
write the effective potential as 



v eff = v - in£tj7 2 , 



(4) 



where the the second term is the centrifugal contribution 
and the first term is the total gravitational potential of 
dark matter, stars, and gas evaluated in the plane of the 
disk z = 0: 



V = V q (vj) + V* + (w,p) + V g (w,vj,t). 



(5) 



The axisymmetric part Vo (137) arises from the mass distri- 
bution of all three components (dark matter, stars, and 
interstellar gas). It yields an angular speed fl(vj) and an 
associated epicyclic frequency k(vj) defined by the gra- 
dient of the specific angular momentum, an expression 
that is also sometimes called the Rayleigh discriminant: 



K 2 = 



1 d 

vj 3 dvj 



(6) 



The role of f2 and k are well known in galactic dynamics, 
and their ratio Q/k is one of the fundamental dimension- 
less parameters in the current theory. 

The quantity V* is the (specified) spiral potential pro- 
vided by the disk stars: 



V* 



-A(vj) cos [imp — , 



(7) 



with A{vj) and $(zu) being the amplitude and radial 
phase of the spiral gravitational potential, both of which 
are regarded here as given functions of galactocentric 
radius vj from stellar density-wave theory. In the con- 
vention of density wave theory, the radial wavenumber 
fc ro = $'(vj) is negative for trailing spiral waves. The 
asymptotic (or WKBJ) approximation of spiral-density 
wave theory assumes a small tilt angle i of the spiral 
arms, i.e., that tani = m/\k^\vj is small compared to 
unity. To justify the use of linear theory of a sinusoidal 
shape factor for the stellar spiral, the radial forcing am- 
plitude of the stellar spiral arms, |fc ro |A, should be a small 



k j sin! 



(9) 



which is not a small parameter because the large factor 
to/ sin i compensates for the small factor F. A physical 
way of stating the same conclusion is that the spiral grav- 
itational field only needs to produce radial velocities com- 
parable to the turbulent or Alfven speeds to have large ef- 
fects (e.g., shock waves) in the interstellar medium. The 
turbulent or Alfven speeds are much smaller than the 
rotational velocities in giant spirals. 

In this paper, we wish to study not only the effects 
of stellar forcing, but the enhancements produced by 
gaseous self- gravity. In 3-D, the gaseous component of 
the gravitational potential, V s (m, ip, z, t) is related to the 
gas surface density, S, by the Poisson equation for a 
razor-thin disk: 



vj dvj 



1 <9 2 K d 2 V„- 



dV, 



dvj I vj 2 dm 2 dz 1 



4irGY,(vj,ip,t)6(z). 

(10) 

In equations ([2]) and (j3)), the radial and tangential com- 
ponents of the Lorentz force per unit mass of the con- 
ducting fluid, J-",^ and F Vl are given by 



1-kYj vj 
27rE vj 



d(wB v ) dB^ 



dz 



dp 



d{wB^ 
dvj 



dB^ 



dp 



(11) 



(12) 



where z <C o is the equivalent half-height of 
the gaseous disk over which the matter is realisti- 
cally distributed. In this paper, we implic it ly as - 
sume that zq is a constant, but iPiddingtonl (|1973[ ) 
pointed out that this state of affairs would lead to 
an enhancement of synchrotron radiation behind spi- 
ral arms t hat is larger than was subsequently ob - 
served (e.g..lMathewson. van der Kruit 1 fc Brouwlll972T) . 



iMouschovias. Shu, fc Woodward! ( 1974 ) proposed that 
magnetic buckling o f the field an d its subsequent infla- 
tion by cosmic rays ((Parker 1969) could solve this diffi- 
culty. In the current analysis, we ignore this complica- 
tion as well as the role of cosmic rays, but we warn that 
more accurate feathering analyses will need modification 
when the feather spacing becomes comparable to the disk 
thickness. 

On the large scales of interest to the problem, the in- 
terstellar magnetic field can be assumed to satisfy the 
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condition of field freezing for a planar magnetic field: 



dB„ i a 

— 5T~ H sr \ B ^ u v - B f u ^> = °! 



as. 







9t C?tx7 



— (B^Up - B v u m ) = 0. 



(13) 
(14) 



Notice that zu times the first equation followed by partial 
differentiation by w added to the partial differentiation 
of the second equation by tp implies that the constraint 
of no magnetic monopoles, 



1 d {wB^) 1 dB^ 



dz 



w dip 



0. 



(15) 



holds for all time if it is satisfied initially. 

Finally, to close our set of basic equations, we model 
the tu rbulent gas pressure with a logatropic equation of 
state (jLizano fc Shul[l98l : 



n g = E u t 2 In 



(16) 



where i>to is a characteristic turbulent speed. The square 
of the signal speed associated with this equation of state 
is given by 



•'to 



(17) 



which mimics the observed tendency for dense interstel- 
lar gas (e.g., molecular cloud complexes) to have lower 
turbulent speeds than rarified interstellar gas (e.g., H I 
clouds). The derivative of n g being positive and decreas- 
ing with increasing £ are more important properties of 
the logarithmic law than the formal feature of having a 
negative turbulent pressure in the regions of low surface 
density, because the formal pressure can contain an ar- 
bitrary addictive constant without having any physical 
effects on the analysis. 

2.2. Axisymmetric State 

We identify the axisymmetric quantities as the zeroth 
order reference state, and denote them by the subscript 
0. With only circular velocities in the corotating frame, 
u v = vj [D,(w) — O p ], and toroidal magnetic fields, B v = 
B^o, that depend on w, the equation for radial force 
balance becomes 



wVL 2 



dv Q i dn 

dm Sg dzu 



z Q B 



27tSq vo dvo 



(ujB v0 ). (18) 



2.3. TASS State 



In the corotating frame, the TASS state is also time- 
steady, so the field-freezing equation ([IB"]) can be satisfied, 
just as in the axisymmetric state, by assuming that the 
magnetic field is parallel to the vector velocity, which we 
can write in the form: 

B m {zu, tp) = &E(n7, (p)u^(ru, tp); (19) 
B v (m, tp) = 6E(nj, ip)u v (m, tp), (20) 

where the scalar factor of proportionally b is chosen 
to be a constant in order to satisfy the condition of 
zero monopoles when the equation of continuity for 
the gas also holds (i.e., B and Su both have zero 



two-dimensional divergence). Because the fluid veloc- 
ity is mostly circular even in the TASS flow, the <p- 
component of the magnetic field is much larger than its 
zu-component, except near corotation where O(ro) = Sl p . 
Far from corotation, if we suppose the asymptotic ap- 
proximation that the TASS flow produces radial vari- 
ations that are large compared to tangential variations 
(or obtained by dividing perturbational quantities by w), 
we may approximate the above expressions by 



Zq 
~2tt£ 



B v d w (B„) ~ -v^-^ Q-) . (21. 



and, 



^^W)^ ^-^-j, (22) 

where we define the square of the unperturbed Alfven 
speed as 



4tt£ /2V 



and we have ignored the spatial variation of So (axisym- 
metric part) in comparison with those of E. 

In giant spiral galaxies, the squares of the character- 
istic turbulent and Alfven speeds, v% and v\ are small 
compared to the square of the flow velocity on the large 
scale, e.g., va 2 Q 2 . In these circumstances, the second and 
third terms on the right- hand- side of equation (|18|) are 
small in comparison to the first, and the rotation speed 
wVl{vj) depends mostly on the gravitational potential 
Vo(nj) of the axisymmetric distribution of dark plus or- 
dinary matter, which we shall henceforth assume to be 
fixed. 

The adoption of the logarithmic equation of state al- 
lows us to write S _1 Vn g = V% g , where H g is the spe- 
cific enthalpy of the turbulent gas: 



Kg = -%) I — 



So 
S 



(23) 



Similarly, we may identify the "specific enthalpy" associ- 
ated with the dominant part of the magnetic "pressure" 
II m = ( V 2 /2)E- 1 S 2 : 



(24) 



Hm — vao \ 



NONLINEAR PERTURBATION 

We now return to the rest of our equations and assume 
that the actual situation is a combination of an axisym- 
metric, time-independent state plus a nonlinear TASS 
response and further feathering perturbations. The sub- 
script 1 in this section will refer to both the TASS re- 
sponse and the feathering instability, but in the Appen- 
dices we shall apply it only to the feathering perturba- 
tions and include the TASS response with the axisym- 
metric state (as unscripted variables when it will cause 
no confusion). In ensuing sections, we avoid confusion 
by attaching a ~ when we mean the perturbations due 
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to the feathering instability alone. Thus, 
S = £0(537) + Ei(t<7, <p, t), 
H = Ho(w) + Wi(t*7, ip, t), 
u-co = u^ {w) + u\{m, ip, t), 
u v = ru(Q — ftp) + vi(zu, (p, rj), 

Ft* = f0(m) + /U7l(ro,<^,t), 

(25) 

where % is the enthalpy associated with the gas pressure 
and magnetic pressure. Note that only v\ is small com- 
pared to its zeroth-order counterpart, with even the last 
approximation breaking down near corotation. Consis- 
tent with the approximation that the pressure gradients 
of the gas turbulent motions and magnetic fields con- 
tribute little to the axisymmetric force balance, we at- 
tribute their influence to the perturbations marked out 
by the subscript 1. Without linearization (because the 
TASS response is highly nonlinear), the substitution of 
equations (|25|) yields the set: 



du 1 



Ui 



dw 



o - a, + — 



dui 

dtp 



w 



-2n Vl --(H gl 



Ha 



(26) 



be constant, then $(tu) = — (m/tanz) la{w/mo), which 
corresponds to the case when the stellar spiral arms are 
fitted by logarithmic spirals. 

We now introduce the orthogonal spiral coordinates 77 
and £ used by SMR with e v x = e z and 



dr\ = — $' ' (vj)dvj + mdip — m ( h dip 

tani w 



(30) 



d£ = 



tani 



m dw 

<I>'(n7) W 2 



dtp 



dw 



w tani 



-d<p 



(31) 



Note that if we move radially outward at fixed w, £ will 
increase by 27r cot i for the same increase in w that re- 
sults in an increase of 2tt for 77. When we transform 
from (w, ip) to (rj, £) and draw rectangular boxes in (r), £), 
the coordinate system is similar to the one defined in 
iKim fc Ostrikerl ([2002) , but there are two differences. 
(1) We do the calculations in a standard Eulerian man- 
ner, without mixing time and space coordinates as in the 
"shearing sheet" treatment. (2) The ratio of the axes 
are depicted in their correct geometric proportions, de- 
termined by the spiral pitch angle i (see Fig. 2). 



dt dw V p w J dp 



U1V1 



w 



2n 



i_d_ 

137 dip 



(H gl +U) 
dH m i 



w(Q — f2 p ) + vi dw 



(27) 



In the above, U = V* + V g is the perturbation of the grav- 
itational potential beyond the axisymmetric state. Con- 
sistent with the approximation made above of ignoring 
the radius of curvature, the continuity equation becomes 



ggi 
dt 



_d_ 

dw 



[(E + S 1 )«i] + ^-[Ei(n-fip)]=0. (28) 



2.4. Spiral Coordinates and Asymptotic Approximation 

We follow iRoberti (|1969f ) and SMR in introducing the 
local orthogonal coordinates (rj, f ) in the plane of the disk 
galaxy, where curves of £=constant and 7y=constant de- 
fine, respectively, the directions perpendicular and par- 
allel to a background stellar density waves with a locus 
of local gravitational potential minimum whenever the 
spiral phrase, 



77(177, ip) = mip — $(177) 



(29) 



that enters in equation ([7]) equals zero or integer multi- 
ples of 2tt. But the function r\(w, ip) can increase by 2tt 
either by ip increasing by 7r (for m — 2) with w fixed (go- 
ing around halfway the galaxy in a circle), or by $(777) 
increasing by 2ir with ip fixed (going out radially until 
the next spiral arm). We wish the solution to look the 
same in either case to lowest asymptotic order. In a local 
treatment, where we approximate the inclination angle i 
of the spiral arms with respect to the circular direction to 



spiral arm 




pitch 
angle i 



galactocentric circle 



Figure 1. Spiral coordinates (r), £) are defined in the direction 
parallel and perpendicular to the spiral arm, locally at a galacto- 
centric radius vj. 



The two coordinate systems are related by the follow- 
ing metric: 



ds z = dw z + w 2 dip 2 



T 2 a\r, 2 i. 



'-(dr, 2 +de), (32) 



which corresponds to a local rotation through an angle i 
and rescaling of lengths by a common factor of ttj sin i/m. 



2.5. Non-dimensionalization 
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Figure 2. A diagram showing the local rectangular box for the 
spiral coordinates (77, £). The spiral arms are indicated by the bold 
vertical lines on two sides. The perpendicular distance between two 
arms is given by harm = L v = 27rro sin i/m, and the dash diagonal 
is the line constant galactocentric radius. Since the solution is 
doubly-periodic in both the r\ and £ directions, coordinates (0, 0) 
and (L V ,L^) represent the same location, and L = cot i = L^/L v . 



After the transformation into the spiral coordinate 
(?7,£), the equations of motion become 



m 



u v i 
dt 


2flu (1 j 


du v i 


d 


' dt 


dr] 



du v i 
dr] 



— (H gl +H ml +U), (33) 



wsini ( dun k 2 \ sdun 



+ (u (Q + 



dt 20 



The corresponding continuity equation reads 

— rjrW So + Sl)Ko + tv)] 

+ ^[(E + S 1 )Ko + «ei)] = 0. (35) 

To write the equations in dimensionless form, we start 
by picking relevant velocity scales for the problem: 

u = u, ql /(2UV) 1/2 and v = u^/V, (36) 
where we have followed SMR by defining 



U 



and V 



20™ 



(37) 



such that the Coriolis terms become v and — u for the 77- 
and f momentum equations, respectively. Similarly, we 
define x 2 = v 2 /2UV and x\ = v\ /2UV for the square 
of turbulent and Alfven's speeds. For the record, we can 



rewrite the enthalpies into the dimensionless form: 

h s = n gl /(2UV) = -x t0 /(l + a), (38) 
h ia = H ml /{2UV) = x A0 {l + a), (39) 

where a = £i/£o is the relative gas surface density be- 
tween the perturbed and axisymmetric states. We also 
rewrite the perturbed gravitational potential U = V* + V g 

as 

U = [-{uj 2 n 2 )Fcosr 1 + (27ra7GS o )0], (40) 

where <f> is the perturbed self-gravitational potential of 
the gas in units of (zu sin i/m)2irGT, . Finally, we mea- 
sure time in units of inverse epicyclic frequency: 



dr = ndt. 



(41) 



We introduce now the following additional dimension- 
less parameters: 



v = -u r]0 /{2UV) 1 ' 2 

/ 



m(£l p — 0)//c, 



(-)• 

\smi J 



(rosin i/m) 2-7rGSo 27rmGSo 

(X = = . 

2UV tdk 2 s'mi 



(42) 
(43) 

(44) 



where / is the afore-mentioned true dimensionless mea- 
sure of the nonlinearity of the stellar forcing, and a is a 
similar dimensionless measure of the strength of the self- 
gravity of the gas. Although 27rGEo may be regarded as a 
small correction to the axisymmetric gravitational field of 
the galaxy, wVt 2 , nonlinear compressions behind galactic 
shocks make gas self-gravity a fierce contractional com- 
petitor to the vortical spinup represented by zuk 2 sin i 
when a is an order unity parameter. 

The continuity and momentum equations now take the 
dimensionless form: 



da 



Y^ 1 



d 
du 

d~T +{ 



a){-v + u)] 

(1 + ff) (-t£7 + ^) 



= 0, 



=v — 



(1- 
dv 

d~r + { - 



20 

u + — 

K 



drj 

Xto da dcj) 



du / v k \ du 
dv + V tan7 + 20 V ~d£ 



dv / V K \ Ov 

'&n + V tan7 + 20 V ~di 



drj 

xto 



tani 
da 



dv 
d£ 



(1 + a) 2 dt dt 



(45) 



(46) 



(47) 



where f v and /j are the components of the dimensionless 
Lorentz force per unit mass in the directions perpendic- 
ular and along the spiral arm, respectively: 



(2uvy/ 2 '' 



_£_ 

V ' 



(48) 



For the closure of the equations, we also need the solution 
of the Poisson's equation for the self-gravity of the gas 
and the equation of field freezing for the magnetic field. 
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3. ONE-DIMENSIONAL SPIRAL SHOCK 

In this section we revisit the steady, 1-D, TASS solu- 
tion, adding in the consideration of the effects of mag- 
netic field (see also Roberts fc Yuan 1970) an d self- 
gravity of the gas (see also iLubow et al.| 119861 ). The 
TASS state, denoted by hats, gives the background flow 
of the feathering problem. All hatted quantities depend 
only on 77, in the form q = q(r]). 

The Lorentz force per unit mass now reads, 



and, 



4 da 



- 2Q ( tani \ da 

h = — Xao tt^ t- 

k \ 1 + a J at] 



(49) 
(50) 



Note that the Lorentz acceleration in t wo dire ctions differ 
in scale by an extra factor of 20,/ k (= y/21/V /V) because 
of the difference in defining the dimensionless u and v. 
Nevertheless, is smaller than f n by the factor tani, 
and it can be dropped asymptotically in the dynamical 
equation for uj. By also dropping the derivatives in r 
and £, the governing equations for velocity in r\ and £ 
directions now read as follows: 

du , ^.v — add/dri—fwxri . 

-r = (-v + u ) — 7 — — — - — ' ( 51 ) 
arj \—v + u) z — x 



du 
dr) 



(52) 



where x = xto(l + 5") + ^Ao(l + &) is the square of 
the effective signal speed, and the parameters / and a 
measure the relative strength of the stellar and gaseous 
perturbation gravitational fields. We make use of the 
dimensionless mass flux as a conserved quantity along the 
flow by integrating the continuity equation and putting 
it into the form: 



(1 + &)(-!/ + U) 



(53) 



The self-gravity term can be obtained from a given sur- 
face density a(rj) by solving the Poisson equation under 
the WKBJ approximation. The derivation is standard 
and given in Appendix The solution can be expressed 
in Fourier series form once a {if) has been found by inte- 
grating the ODEs (|51l and 1521 for u and v, coupled with 
equation 1531 for a): 

00 

a = + cos ( nr l) + S n sin (nrj)] , (54) 



n=l 



and 



c!4> 
dr) 



[-S n cos (nrj) + C n sin (nrj)] , (55) 



where C n and S n are the n-th Fourier components for 
even and odd solutions, respectively. On the other hand, 
the integration of equation (|51[) requires knowledge of 
4>, so iteration (with a relaxation parameter) is required 
to find numerically a completely self-consistent solution. 
Apart from the added iteration and convergence steps 
for the self-gravity when a is nonzero, and a different 
expression for x, the equations (|51I52[) have the same 



form as the set studied by SMR, and they can be solved 
by using the same shooting method with the matching 
of upstream and downstream flows satisfying the shock 
jump conditions discussed below. 

3.1. Magnetosonic Point and Shock Jump Conditions 

The spiral shock solution is periodic in the r) direction 
in the sense that when the flow passes through the shock 
front, it will accelerate from the submagnetosonic speed 
to supermagnetosonic speed, and eventually reach an- 
other shock at the next spiral arm. Thus, the region be- 
tween two consecutive shocks is transmagnetosonic and 
has a magnetosonic point location (77 = i] mp ), where the 
speed of the flow equals the local speed of magnetosound. 
Solutions with multiple magnetosonic points and shocks 
are also possible, but th eir study is beyond the s cope of 
this paper (see SMR and lChakrabarti et al.ll2003l for dis- 
cussions of the role of ultraharmonic resonances for pro- 
ducing spiral branches and their possible relationship to 
flocculence when overlapping resonance leads to chaotic 
nonlinear behavior). The magnetosonic point is located 
where the following condition is satisfied, 



{-v + uf 



0. 



(56) 



which is also an apparent singular point of the equation 
(|5"Tj) . By substituting equation (|53")) and the equation of 
state, we get 

{-v + iif + — (-is + ii) 2 + x A0 v = 0, (57) 

V 

which is a cubic equation that gives only one positive 
value of ( — v + u) = {—v + ii mp ) algebraically if v is 
negative. A smooth solution across the magnetosonic 
point can be found by requiring both the numerator and 
denominator to be zero in the equation (j51[) . Therefore, 
the derivatives of u and v at the magnetosonic point can 
be evaluated as: 



du 
dr] 



[-Ump/y, 



mp 



Imp 



/ cos rj 



mpj 



11/2 



[2 + xtov 1 /y mp - XAov/yi 



1 1/2 



imp] 



dv 
dr] 



(58) 
(59) 



in p 



^mp 



where <f>" is the second ry-derivative of (j>, and we define 
y mp = — v + u mp . Note that the value of the deriva- 
tives can be evaluated by solving u mp in advance from 
equation (I5T1) with the background parameters given. In 
our implementation of the shooting method, we start the 
integration from the neighboring points of the magne- 
tosonic transition to the shock front in both supermagne- 
tosonic and submagnetosonic directions separately. The 
"initial" values of u and v at these points are given by 
the following Taylor's series: 



du 

U = Ump + ^ 



dr\ 



(V - Vmp) + 



nip 



/ sin n mp 



dv 
dr] 



(60) 



(V-Vmp) + 



The derivatives of self-gravitational potential (f> are ob- 
tained from the solution of previous step. Since d4>/dn is 



Lee and Shu 



generally a smooth and continuous function, the value of 
0"|mp may be expressed in a finite difference form with 
little numerical error. 

3.2. Matching Conditions 

The physical problem is constrained by the fact that 
the downstream and upstream flows for a periodic solu- 
tion must match the values that allow a shock jump con- 
ditions to connect the supermagnetosonic and submag- 
netosonic collision. These jump conditions are obtained 
by requiring the sum of gas (turbulent) and magnetic 
pressures and momentum fluxes to be continuous across 
the shock: 

2 



(f + a)(-v + u)u + x w ln(l + a) + ^(1 + a) 2 



and, 



[(l + a)(-u + u)v]l 



to vanish separately. We can identify the constant mass 
flux, (f + <r)(— v + u) = —v from the continuity equation 
above. Thus, we may obtain the corresponding post- 
shock (or pre-shock) values of u and v for a given pair of 
values on the other side of the shock. In practice, we cal- 
culate the corresponding post-shock (submagnetosonic) 
values by using the pre-shock (supermagnetosonic) val- 
ues, as if they were to satisfy the shock jump conditions. 
We postpone the discussion of the numerical results until 

m 

4. FEATHERING ANALYSIS 

In the context of the feathering phenomenon, we need 
to consider variations in 2-D and time. Because the ref- 
erence TASS state is independent of £ and t, we may de- 
scribe, in a linearized treatment, the additional feather- 
ing variations as oscillatory disturbances in £ and r. Such 
a treatment constitutes a standard linear-stability anal- 
ysis and should be contrasted with prior treatments that 
supposed feathering to be a shearing, time-dependent, 
phenomenon that is imposed by the mathematics of a 
transformation that is useful only near corotation. In 
the current paper, we deliberately stay away from coro- 
tation. A complication that does appear is the oscilla- 
tions introduced by a wiggling shock front, which leads to 
perturbed jump conditions that further affects the down- 
stream flow. In any case, instead of solving a set of PDEs 
as in the numerical experiments (which do not need lin- 
earization), we obtain a set of ODEs for the feathering 
perturbation. Imposing double-periodicity, for given m 
(= 2 in the usual TASS picture) and I (the number of 
feathers strung out along the arms in the ^-direction 
per spiral arm box), the (complex) oscillation frequency 
u>r + iuji in time becomes an eigenvalue of the overall 
problem. 

4.1. Perturbational Equations 

As the feathering perturbations are time-dependent 
and vary along both r\ and £ directions spatially, we de- 
fine the variables as follows: 

u = u(rj) +u(r),£,t), 

v = v{rj) +v(i],€,t), 

a = a(r,) + ~a(r 1 ,^t), (61) 



where the hat states are the background TASS flow, and 
the tilde states are perturbations assumed to be small 
compared to the background. The perturbational mag- 
netic field is time-dependent, and we no longer assume 
its direction is parallel to the flow as was assumed for 
the TASS background. We will derive the perturbation 
induction equation in Appendix [B]), where we show it 
corresponds simply to the conservation relation for the 
magnetic flux function A (^-component of the vector po- 
tential for the magnetic field). For here, we simply record 
the resulting linearized perturbation fluid equations for 
the tilde quantities: 



da 

d~r + ^ 



L-, — h Vt 



= -(1 



da da 
an d,t] 

„ dii du _ 

dn dn 



da 



K dv 
'2W (62) 



du 



du 



dii 



+ (—v + u)— — h u- — h i>T 



XtO 



{l + af 



dq 
da 
drj 



dr) 

da 



du 
"dl 



(63) 



2xt 



dv , „.dv 
d-T + { -» + U) d-n 



2ft 

K 



dn (1- 

_ dv 
dn 

da 



vt 



dv 



(1 + a f d£ di 



+ ft, (64) 



where wt is the total ^-component of the background 
fluid velocity (i.e., axisymmetric plus TASS) in the frame 
that corotates with the stellar spiral. The perturbed 
Lorentz acceleration, and are given by 



d 



d 2 



fn = xao TT~^ + M 



drj 2 d^ 



and, 



2tt 



-XAO 



2n 

£A0 



tani 
1 + a 



drj 2 



dAi 



1+a dr] 



(65) 



de 



At 



l + a di 



(66) 



respectively. In the above, A\ is the dimensionless 
z-component of the perturbational magnetic potential. 
The perturbational evolutionary equation for it reads 
(see Appendix B): 



dA x dA t 
~dV + { ~ V + u) ~dn~ 



vt 



dAi 



= (l + a)u- 



(67) 



4.2. Perturbed shock jump conditions 

The perturbed shock jump conditions can be obtained 
by linearizing the jump conditions in the frame of per- 
turbed shock front. The shock front is displaced and no 
longer parallel to the spiral arm. In Figure |31 we show 
the configuration of the perturbation. Inside corotation 
radius, the flow is entering the shock from the left in the 
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pre-shock 



post-shock 



A=2ttL// 



Figure 3. Corrugation of the shock front. The perturbation, e, 
along the shock front is characterized by the wavelength A = 2irL/l 
and is assumed to have small amplitude. The vertical dotted line 
is equilibrium shock front at r\ = rjg^. The unit normal and unit 
tangent are denoted by h and t, respectively. 



frame of the shock. Hence, we must obtain the normal 
direction of the shock and the shock velocity in the cur- 
rent frame. The position of the perturbed shock is now 
given by, 

r) = Vsh + e(t,t), (68) 

where rj s h is the unperturbed position and e is a small 
number. We define, 



s = <q- rj sh - e(£,i), 



(69) 



to be the displacement from the moving shock front. 
Thus, the locus of the shock front is s = and the unit 
normal n is given by, 



1 



|V*| 
1 



ds 
drj 



ds 



(70) 



where |Vs| 



1 + 0(e 2 ) is the magni- 



L + (9e/90 2 ] 1/2 
tude of the gradient normal. The velocity of the shock, 
D, which is normal to the shock front, can be found 
by considering a normal displacement Ar of the locus 
at time At, i.e., |Vs|Ar + At(ds/dt) = 0, and thus, 
D = hAr/At = (de/dt)h in the first order of e. 

There are five shock jump conditions in the problem 
(see ?, eqs 25.16-18,20,21). The linearized perturbation 
in the moving shock frame reads, 



[p5u±_ + u±6p] 1 = 0, 
u 2 ± 5p + 2pu±5u± + 



_ Bj_ 

47T 

[SBx]l 
[B ± 5u 



SB n 



pu±5u\\ 



47T 



-SB 



dP 
~dp~ 

puuSu 

1? = Q, 



6p- 



47r' 



'-SBn 



= 0, 



0, 



B«S 



\\ou± 



l SB ± -u ± SB [l ] 2 1 = 0, 



(71) 
(72) 

(73) 
(74) 
(75) 



where the variables with and without 5 are in first and 
zcroth order of e, respectively. The jump conditions are 
given in terms of variables parallel and perpendicular to 
the shock front locus, s = 0. Since the perturbation on 
the shock front is in first order of e, the non-5 variables 
are simply their background counterpart. For the vari- 
ables that are first order in e, we include both Taylor's ex- 
pansion at the perturbed shock front and the geometrical 
projection (i.e., in Lagrangian sense). Thus, by making 
use of the unit normal in equation (|70|) . and the corre- 
sponding unit tangent, we express the dimensional gas 
surface density, flow velocities and magnetic fields into 
the following, 



p ~ p + p + edp/dr], 



Uj_ 



u 



v 



v -r edurj/dr] — D v — u^d^e, 



uii — U{: + ii^ -f- edu^/drj + u v d^e, 
B ± ~B ri + B^ + edBr,/dr) - B&e, 
B\\ ~ Bj: + B{ + edB^/dri + B^e, 



(76) 
(77) 
(78) 
(79) 
(80) 

where we evaluate the variables at 77 = rj s h, and we de- 
fine D v and as the r\ and £ components of the shock 
velocity, respectively. Note that D v = (de/dt)/\\7s\ 2 ~ 
(de/dt), and we exclude — 0(e 2 ). To make the equa- 
tions dimensionless, we measure the surface density, ve- 
locities and magnetic fields in term of £0, V2UV and 
Bfo, respectively. We now write the boundary conditions 
(|7l1- 175|) a more compact form, 



Q(1)V(1) = Q(2)V(2), 



where Q = Q(r/) is a 5 x 5 matrix is given 
by the coefficients in the equations (|76H80[) in 
terms of the background variables and V = 
[5(7, 5u, (k/2Q)8v, 5B±/B^ , 5B\\/B^] t is a column vec- 
tor of the 8 variables. 

4.3. Stability analysis 

As the governing equations for the perturbed variables 
do not have explicit dependence in r and £, we can sim- 
plify the equations by assuming that the tilde variables 
have e lUT - ll t/ L dependences. The perturbed shock front 
must have the same sinusoidal dependence in time and 
space (see Figure [3]) . The instability condition follows 
when Im(cj) = u)\ < 0. To treat the perturbational self- 
gravity, we use a simplified solution to the Poisson equa- 
tion obtained in Appendix [AJ 



L 



(82) 



appropriate to the assumption (motivated by the obser- 
vations) that the important feathering corresponds to 
large I (many individual feathers in the £ direction for 
the equivalent spiral box where we only have 1 TASS 
arm becoming another in the 77 direction). 
After rearranging the terms, we get 

=(-u' - ilo t )(7 Ui i - a'u u j + ^7^(1 + °)v u ,h (83) 
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XtO 



(l + a)* 
+ (-u + u) 



da^ 



drj 



dr\ 1 + <7 drj 



dA 



lu.l 



dr) 



2x w a' 



(1 + <r) 3 



Pa, I + (—u' - iuir)u u ,l + V u ,l 



(84) 



•''.\n ( j ) *h^..i. 



^dv u ,i 20 

{-V + U)— 1 XAO 

drj k 



tanz \ dA' lu , 



20 

K 



(1 + a)u u> l - iUTV u ,l 



a J drj 



Oi~. 



20 



-^AO 



tani / I 



1 + <t \L 



il 



1 + <J \L 



Aiu 



(85) 



where we define ujt = — (1/L)vt to be the dimension- 
less Doppler-shifted frequency in the moving frame of the 
background flow along the spiral arm. The transformed 
induction equation reads, 



{-u + u)A' lu)jl 
= - iujTAiuj + (l + a)u u ,i 



(86) 



(— 

V20 



tani v u ,i 



Similarly, by taking the Fourier transform and matching 
the first order terms in the perturbational jump condi- 
tions (|76H80p , we can express the 6 terms in the in terms 
of the tilde (Eulerian) variables: 



5a 
5u 

Sv = 



5B ± /B i0 
SB l{ /B C0 



a + eda/drj, 

u + edu/dr) — iuiTe, 

, A , , 20 il „ 

v + edv/dri -ue, 

K L 

il ~ il , 
_ Al + j(l + a )e, 

-A[ + eda/dr), 



(87) 



where we drop the subscripts (w, I) for clarity, and we 
have used dB v /drj oc d(T,u v )/dri = 0, 



El 



e dBs _ e d£ 
B iQ drj S drj 

= (l + a) and 



da 
drj ' 



B 



(0) 



B„ 



tani, 



B^o B^o B^ B^ 

and we neglect the term with d^e(B v / B^q) = O(etani). 

4.4. Method of solution 

We try to solve the four perturbed equations (|83ti86[) 
as a set of ODEs. Since the Lorentz force contains a sec- 
ond derivative of the scalar magnetic potential, it might 
seem that we require one more differential equation for 
its derivative, dA\/dr) = A[. However, the induction 



equation (|86l) is in fact an algebraic equation for u, v, A\ 
and A'i, and does not involve dA\jdr). The physical rea- 
son is that field freezing implies that the tilde magnetic 
potential must yield a magnetic field structure that cor- 
responds to the TASS magnetic field stretched in time by 
the motion of the electrically conducting matter in the 
feathering instability (see Appendix B). In other words, 
the field freezing equation (|86[) must hold in space simul- 
taneously with the other differential equations, and so, 
the system is a set of Differential Algebraic Equations 
(DAEs), where dA\jdr\ is obtained by differentiation af- 
ter we make the set self-consistent by treating uj not as an 
arbitrary constant, but as a (complex) eigenvalue of the 
time-dependent transformation of the TASS state to one 
that contains (exponentially growing) feathering pertur- 
bations. Fortunately, the problem so posed can be solved 
by following the procedure discussed below. 

First, we should reduce the order of the perturbed 
equations such that it solves for the tilde variables 
(a, u, v, Ai) (and their first derivatives on the RHS) only. 

The solution of A[ can be obtained algebraically from the 
induction equation once we have solved for other vari- 
ables. The elimination of A" may then be done by nu- 
merically differentiation of the induction equation, i.e., 

, ^du u< i k .dvu.i , , , -s^'iw.i 

- (1 +a) ^r + 20 tani ^ + { - v+u) ^trr 

j) (^j) (1 + &)Av*,i - {& + tvr)AL,i- 



Second, we should reduce the number of perturbational 
jump conditions to four, as we have four differential 
equations after the reduction above. In fact, the fifth 
jump condition, equation (|75[) . is satisfied automatically 
by the "algebraic" induction equation. In other words, if 
the complex frequency w has the correct eigenvalue, we 
can impose both double-periodicity and match all the 
requisite jump conditions, with the equation (|75|) being 
consistent with the induction equation with which we use 
to calculate A[. We can eliminate A[ in the connection 
conditions by using the equation (|86[) , which is valid for 
both sides of the shock separately. Therefore, the di- 
mension of the coefficient matrix Q in equation (181 1) is 
reduced to 4 x 4. 

After we "separate" the induction equation from the 
differential equations in the above manner, we may solve 
the system as a Two-Point Boundary Value Problem with 
Eigenvalues. St andard methods of attack exist in the 
literatures, e.g., lAscher et al.l ([1995). One difficulty of 
using the publicly available numerical packages is that 
they do not treat the embedded jump conditions present 
in our system. To solve this problem, we can artificially 
modify our jump conditions in the following form: 



Q(l)V(l) = d and Q(2)V(2) = C 2 , 



(89) 



where the two vectors Ci and C2 are varied until they are 
equal to each other. The procedure leads to the ability 
to use standard packages at the numerical expense of 
solving four more equations. There are also standard 
methods to solve a system with complex variables as in 
our problem, and we do not discuss the numerical issues 



Feathering Instability of Spiral Arms. I. Formulation of the Problem 



11 



Table 1 

Typical galactic parameters for the feathering example of this 
paper 



Rotation curve 
Galactocentric radius a 

Farm 

Turbulent speed 
Alfven's speed 
Magnetic field b ' c 
Mean gas surface density c 
Inclination of stellar spiral arm 
Pattern speed c 



Sofue et al. (19991 
5.0 kpc 
3.8 kpc 
13.4 km/s 
13.4 km/s 
14.8 fiG 
38.6 M /pc 2 
14° 

23.4 km/s/kpc 



a the distance from the modelling region to the galactic 
center 

b value of B,po for scale height zq = 200pc 
c for a = 0.35 
d adopted from ' 



adopted from 



Kendall ct al. (2008) 
Westpfahll iWm) 



Table 2 

Typical set of local parameters (vd = 5.0kpc) for the feathering 
example of this paper. 





0.666 


tan i 


0.249 


F 


11.5% 


a 


0.35 


V 


-0.666 




0.1 


ZtO 


0.1 



further here. 



5. NUMERICAL RESULTS 



In this section we present resu lts of the calculatio n with 
input parameters based on the iSofue et "all (|T999) rota- 
tion curve of the M81 Galaxy. Table Q] lists the relevant 
numerical values used for the feathering calculation, and 
Table [2] gives the seven dimensionless parameters needed 
for solving the governing equations in that calculation. 
The listed mean gas surface density E , gas turbulent ve- 
locity i>to, an d interstellar magnetic field are all too high 
for a Sb spiral galaxy of luminosity class II like M81, but 
we adopt these extreme values to illustrate a point that 
will become apparent in our closing discussion. 

In any case, the mean plasma /3 is /3o = Xto/^Ao = 1 for 
the choice Xto — xaq, and implies that turbulent stresses 
cx Xto/(l + &) dominate in the inter-arm region where 
a < 0) while magnetic stresses oc xao(1 + 0") dominate in 
the arm region where a > 0. We define a mean Toomre's 
Q parameter for the gas accounting for mean turbulent 
motions and magnetic field (with mean effective sound 
speed, ao oc y^Eo) as: 



Q 



Kao 



-(sGtO + £A0 



,1/2 



(90) 



With Xto — xao — 0.1 and a — 0.35, we have Q = 2.55, 
so the gas is stable on average to all axisymmetric self- 
gravitational perturbations. 

5.1. TASS profiles with self-gravity and magnetic field 

The dimensionless parameter a characterizes the 
gaseous self-gravity. In general, it is an order unity quan- 



tity for spiral galaxies of not too early a Hubble type. 
Figure [4] shows the surface density in the perpe ndicular 
direct ion to the arm. The plot is similar to the iRobertsI 
(1969) calculation, except the surface density is no longer 
arbitrarily scalable when we include self-gravity. The 
shock strength increases with increasing a (with all other 
parameters held fixed as in Table 2) because the gaseous 
self-gravity deepens the minimum of the spiral gravita- 
tional potential. More of the support for the total spiral 
gravitational potential coming from the gas also pulls the 
shock front downstream closer to the potential minimum. 
Increases in a also rounds out the density peak. These 
effects were also seen in the 1-D numerical simulations of 
IKim fc OstrlkeH ((20021) . 



150 



100- 




Figure 4. Surface density of the gas for a from 0.1 to 0.35 (cor- 
responding to 11 to 39 Mq/pc 2 for the mean surface density) and 
F = 11.5%. The minimum stellar gravitational potential is located 
at rj = 0. The thick line (a = 0.35) is the background density pro- 
file used for the feathering perturbation example. 



On the other hand, the full-width-at-half-maximum 
of the density profile is approximately constant. If the 
FWHM is taken as a representative value, the width of 
the arm with spiral galactic shocks is about 12% of the 
distance between two arms, or about 480 pc in the cur- 
rent model. 

MAGNETIC FIELD 




Figure 5. Surface density of the gas for plasma j3, /3rj = 1.0, 2.0, 10 
and 20 (corresponding to B v o = 7.9, 5.6, 2.5 and 1.8 fiG, re- 
spectively). The mean surface gas density is set to 11 Mq/pc 2 
(a = 0.1). 
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Co ntrary to previous assertions (cf. iDobbs fc BonnelH 
2006), the magnetic field plays an important role in spi- 
ral galactic sho cks and the resultant f eathering instabil- 
ities (see, e.g.. iKim fc Ostrikerl 120021 ). Figure [5] shows 
some different choices for the magnetization parameter 
xao = Zto/Po, when the turbulent and self-gravity pa- 
rameters are kept fixed at xto = 0.1 and a = 0.1 with all 
other dimensionless parameters held at the values given 
in Table 2. In general, the increase in xao (or decrease in 
/3o) suppresses the compression of gas in the postshock 
region, and lowers the shock strength. Conversely, the 
peak surface density rises very rapidly with increasing 
Po (i- e -i decreasing magnetic field strength), and there 
is no steady solution possible for (3q much larger than 
20; i.e., the spiral arms would go into continued gravita- 
tional collapse with a = 0.1 if the magnetic field is too 
weak. The magnetic field cannot be ignored either for 
the structure of the TASS density pattern or for the de- 
velopment of the feathering instability when self-gravity 
is important. The dependences of the background pro- 
file and the feathering perturbation on the dimensionless 
parameters of the problem will be investigated in Paper 
II. 

5.2. Feathering Perturbation 

Taking the TASS 1-D solution as the background for 
the feathering phenomenon, we solve the perturbed equa- 
tions for each Z-mode and obtain the 2-D solution using 
inverse Fourier transforms. A larger survey of parameter 
space is undertaken in Paper II. Here, we just show a 
typical result for the total surface density and magnetic 
field lines in Figures |5] To obtain sufficient contrast, we 
have arbitrarily scaled the linear perturbations so that 
they are no longer small compared to the background. 
Figure [6] compares the flow solutions with and without 
the feathering perturbation for 1 = 8. For better view- 
ing of the post-shock region, the horizontal axis is not 
really n, but n — ?7 s h. For the background flow on the 
left, the magnetic field reaches peak compression behind 
the shock but become weaker for increasing rj as the ex- 
pansional flow out of the spiral arm pulls apart the the 
frozen-in field lines. With the development of the feath- 
ering perturbation with ( = 8 on the right (i.e., 8 feath- 
ers in a distance, Lg), over-dense regions jut out from 
the spiral arm toward the interarm region downstream. 
In this particular case, the Doppler-shifted frequency is 
wt = —0.113 — 0.174i in unit of k. The non-zero real 
part of wt implies that the pattern of feathers moves 
along the outward £ direction with the passage of time, a 
result also seen in the numerical calculations of the Os- 
triker group. The negative value of the imaginary part 
of ojt implies that this mode is unstable and can be ex- 
pected indeed to grow to nonlinear amplitudes with the 
passage of time. With n = 70.2 km s _1 kpc -1 , the e- 
folding growth time-scale t g is 



to 



1 



1 



KlmujT 0.174k 



80 x 10V- 



(91) 



In Figure [71 we show the velocity field (arrows) super- 
imposed on the surface density profile (colored contours) 
of the feathering perturbation at a single instant of time. 
The perturbed flow follows closely to the background spi- 
ral shock profile because the background circular motion 



is dominant over all other motions. Nevertheless, signif- 
icant convergence toward density peaks and divergences 
from density troughs can still be found along the spi- 
ral arm, especially toward the beginning of the feathers. 
The be havior can be profi t ably c ompared to the zoom-in 
plot of Shetty & Ostrikcr (2006). It is tempting to spec- 
ulate whether the velocity fluctuations in the nonlinear 
development of the instability can be mistaken for turbu- 
lent velocities in insufficiently angularly resolved images 
of spiral galaxies. 
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Figure 6. Comparison between the background flow with a = 
0.35 (left) and the flow with feathering perturbation of the I = 8 
mode (right). The color shows the surface density of the gas in 
linear scale, and the white lines are the magnetic field lines. Since 
L = cot i ~ 4 in the model, the periodicity of £ is approximately 
8tt. 



6. DISCUSSIONS AND CONCLUSION 

It is illuminating to compare the numerical result (|9 1 1) 
with a simple model of 1-D Jeans instability. The Jeans 
instability cannot occur in 1-D equilibrium states if the 
compression occurs isothermally because the increase 
in pressure for ces keeps pa ce with the increase in self- 
gravitation fcf. lSpitze rl ll96 8). but if the (turbulent) equa- 
tion of state is softer than isothermal, as it is for the 
adopted logatropic law of the current paper, then it is 
possible for the self-gravity of condensations parallel to 
the galactic shock to overwhelm the declining resistance 
of the tur b ulent motions. According to the estimate by 
iShu et all (|2007l ). the feathering instability is then basi- 
cally the self-gravitational contraction of over-dense gas 
along post-shock magnetic field lines, which is almost 
aligned along the density-ridge of the TASS pattern. A 
rough estimate of the contraction time is then 



t c =s T 



GUgW 

L 2 



-1/2 



(92) 
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Figure 7. Velocity field of the feathering perturbation (I = 8) 
with the background spiral shock. Without feathering perturba- 
tion, the velocity is expected to have no variation along £. 



where W is the half-width of the TASS spiral arm, which 
we shall take to equal the 0.48 kpc mentioned previously, 
and L — L arm (L/l)/A is the quarter-wavelength of the 
feathering instability (the center to edge distance of over- 
dense regions) and also equals, by coincidence, 0.48 kpc. 
If we replace E ff with the mean surface density in the 
spiral arm, which is a few times denser than the average 
including the interarm region, for I = 8 and a = 0.35, we 
have 



(4.302 x 10~ 6 )(71 x 10 6 )(0.48) 



(0.48) 2 
= 0.031(km/s/kpc) _1 
= 31 x 10 6 yr, 



-1/2 



(93) 



where we have chosen E s to have its half-peak value, 71 
Mq pc -2 , in the feather. 

The rough estimate (l93l) underestimates the accurate 
computation of equation (|9"Tj) by a factor of 2.6, suggest- 
ing that differential expansion and magnetic stresses in 
the postshock region play stabilizing influences in the ac- 
tual feathering phenomenon. Nevertheless, although the 
correct mathematical calculation of the feathering insta- 
bility is complex, the basic mechanism behind its opera- 
tion is roughly quasi 1-D contraction along magnetic field 
lines roughly parallel to the spiral arm while the back- 
ground flow is swept downstream roughly perpendicular 
to the spiral arm. Numerical simulation then demon- 
strates that the full nonlinear development of the insta- 
bility results, not in permanent collapse for most of the 
gas in the over-dense regions, but to a redispersal in be- 
tween the spiral arms as the expanding set of background 
magnetic fields helps to tear apart the dense condensa- 
tions that might in a nonmagnetic context h ave experi- 
enced overall continued gravitational collapse. iShu et all 
(2007) suggest that this is the reason why OB star for- 
mation, as prominently as they seem to delineate spi- 
ral structure and substructure, is actually quite ineffi- 
cient in its operation in the present, relatively strongly 
magnetized, universe of interstellar media. This men- 
tal construct, coupled with the visual display of Figure 
[71 suggests that giant molecular cloud associations are, 
not permanent material entities, but a manifestation of 
the parasitic formation and dissolution of the feathers 
at the crests of the nonlinear density waves that we call 
gaseous spiral arms. With this point of view (which we 
recognize will not be universally accepted), the pattern 
is long-lived; the individual condensations are not. 

This research is part of the PhD Thesis Dissertation 
of WKL in the Physics Department of UCSD. The au- 
thors would also like to acknowledge the support of the 
National Science Council (NSC) of Taiwan for its sup- 
port of the Theoretical Institute for Advanced Research 
in Astrophysics (TIARA) based in Academia Sinica's In- 
stitute of Astronomy and Astrophysics (ASIAA). 



APPENDIX 

WKBJ APPROXIMATION OF SELF-GRAVITY 

The self-gravity of the gas is governed by the Poisson equation in the thin-disk geometry, 

V 2 V g = At:GT,5(z), (Al) 

where V g and £ are the gravitational potential and surface density of the gas, respectively. We use the above equa- 
tion for both quasi-lD spiral shock and the 2D feathering perturbation. In the asymptotic approximation (quasi- 
rectangular) of the spiral coordinates introduced in this paper, we can write the Laplacian in the following form: 



V 2 = 



za sim 



dC, 2 



(A2) 



where C, = z/[(zu sin i/m)] and z is the physical coordinate perpendicular to the plane of the razor-thin disk. Consistent 
with the equation (|4T))) . we define V g = 2ttGY,o(vu smi /m)(j>, so that </> is the dimensionless gaseous self- gravitational 
potential. The dimensionless Poisson equation now reads, 



drj 1 



d( 2 



(A3) 



and is subject to the following boundary conditions: 



2a as the integrability condition across the midplane 



( = 0; 4> when £ — ► ±oo; periodic boundary conditions for rj and £ directions. Note that both <fi and dr,<j) are 
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continuous across the spiral shock, but the term d 2 <f> requires special treatment because of the delta function on the 
right-hand side. If we Fourier transform in 77 and £, we have [— n 2 — (l/L) 2 + d 2 ](f> n .i = 2a n ,i8(£), where n and l/L 
are the corresponding wavenumbers for these directions. By integrating the last expression across C — and requiring 
exponentially decaying solutions in the £ direction for both positive and negative values of £, we get the familiar WKBJ 
result: 

\,,= , (A4) 

'■ + {i/ly 

For a 1-D TASS density profile, we take / = and get <p ny o = — 2a n> o/\n\. Then it is straight forward to obtain the 
equations (|54j) and (|55|) by considering the real and imaginary parts of a n . 

For the 2-D feathering perturbation, we adopt a simplification of the complete treatment. For feathering pertur- 
bations reminiscent of the substructures observed in real spiral galaxies, l/L is appreciably larger than the n values 
needed to reconstruct a reasonable accurate surface density profile of the TASS background state. To be sure, a 
formally infinite number of n's are required if we wish to recover the sharp jump of the background shockfront, but 
this is a feature of the background state and not of the smoother perturbations that we are ascribing to the feathering 
instability. Correcting for the finite thickness of the disk would also lead to smoother relations between the perturbed 
self-gravitational potential and the perturbational surface density. With the assumption that l/L is much larger than 
the n in any of the important Fourier coefficients <f) n> i and cr,, j, we adopt the following approximation: 

&(*7) = -777frSi(»7) and 4>'i(v) = -TTTfT^'fa)- ( A5 ) 

\l/L\ \l/L\ 

In the actual example shown in Figure the validity of the approximation is questionable, as it amounts to the 
assumption that (l/L) 2 = 2 2 is a large number. But the feathering displayed in that figure also has too large a spa cing 
between condensations; better examples will be given in Paper II where the approximation used in equation (|A5j) has 
more justification. 

PERTURBATION ON THE LORENTZ FORCE AND INDUCTION EQUATION 

Because the interstellar magnetic field is free of monopoles, it is derivable as the curl of a vector potential A. For 
a field that lies entirely in a plane, the vector potential can have a single component, A = A(w,ip,t)e z . Under these 
circumstances, 

B = Vx (Ae z ) = -e z x VA (Bl) 



The equation for field freezing can now be written, 



V x 



dA 

—e z - (e z x VA) x u 



= 0. (B2) 



With a proper choice of gauge (namely an initial time-independent state in which B is parallel to u, or VA is 
perpendicular to u), we can "uncurl" the above equation, expand the triple vector product, and derive an evolutionary 
equation for A: 

dA 

— + u-VA = 0. (B3) 

In other words, field freezing in this context is simply the statement of the conservation of A as we follow the motion 
of fluid elements. 

If we write A = Aq + Atass + Ai, and u = Uo + utass + u ij the satisfaction of the condition of field freezing by the 
zeroth order axisymmetric and TASS states implies to linear order that 

dAi 

— + (u + utass) • VAj = - [V (A + Atass)] • ui. (B4) 

The elimination of one spatial derivative in the evolutionary equation for A\ (or A\) explains why there is no equation 
for dA'Jdrj in i fOl 
If we use tilde to denote dimensionless perturbational variables, 

a m A\ 

M = —-5-, (B5 

in sin 1 ij^o 

we then get from equation (|B4|) : 

wsaxidAx dA x dAi rosini 

KT+ U v-^ l " M C _ S^ = - B„U£i) . (B6) 

m at or\ at, m 
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Again, we have used unscripted variables to denote the axisymmetric state (denoted with a zero) plus the TASS 
value (denoted with a hat), and a subscript 1 to denote the dimensioned quantity associated with the feathering 
perturbations (when nondimcnsionalized, these are given a tilde). Using the definition of A\ in equation (|B5|) and 
dividing equation (|B6|) by y/2UV, we obtain the dimensionless induction equation: 



Id Ax , ftAi ( 

+ {-v + u) 



K \ dAi ( K \ _ . 

^j^ = (l + ^-(^tan^, (B7) 



k dt drj \tani 2Q 

where we have used u = u n \ /Vzuv, v = u^i/V, B^/B^o = (1 + a) and B v /B^q = tani. In practice, we also need the 
governing equation of A[ = dAi/dr), which can be obtained by taking the ^-derivative of equation (|67[) : 

1 dA\ , du 7, — v dA\ k dv 8A\ , da _ / k \ dv 

k~W + { - v + u) ^ + ^ + t^i^i + mTrtli = (1 + a) d~ v + T, u - (20 H ej- (B8) 

To compute the linearized Lorentz force, we use the expression 

2z 



4tt£ 



[(V 2 A) VAx + (V 2 Ax) VA] . (B9) 



Again, we have used A as a short hand for A$ + Atass- We also ignore the linearized contribution that comes from 
expanding E = So + Etass + ^1 hi the denominator because we are interested the feathering effect in the postshock 
region, where the perturbation surface density Ei is very small relative to the axisymmetric and TASS contributions. 
More explicitly, then, we have to lowest asymptotic order for small sini: 

, „ m fdA^ OA \ 

V (A + Atass) - — n~e^ + — e ? = -B^e v + B v e^, 

zuswt \or) o£ J 

and, 

, . . . m ( dBf dB„ \ m dBc 

V 2 (A + Atass — -- ^ + ^r)= : — r_ T"^' 

w sin 1 \ orj oc, / ^ sin z 0(77 

where on the right-hand-sides we have used A = Aq + Atass as a shorthand. Therefore, the perturbed Lorentz force 
per unit mass can be written as, 

2za <V 2 A 1 )VA-^-{V 2 A)VA u (BIO) 



4?rE v ' 4vrE 

where we define the perturbed magnetic field, Bi = — e z x VAi. For notational convenience, we write the perturbation 
Lorentz acceleration as coming from two parts: fi = fW + f' 2 ), where 

2z Q 



f (i) = *| (y2 Al) _ s vk 



in 



rosin! 




"iol^ + ^l [e,-^e e ], (Bll) 



and, 




to sin i 1 + (j dr) \ drj 



(B12) 



In the above, we have approximated cos 2 i ~ 1 for small sini. The dimensionless components of the Lorentz force 
(/»)! /<;) can be found by rearranging the terms: 



_ / wsini/m > 
/?? ~ I 21^ I 1? ' 



, a 2 <9 2 \ j 9Aj 
=^ao tt^j + tttot Ai + xao- 



9ry 2 <9£ 2 / 1 + a drj 
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and, 



ft 



■uj sin i/m \ 

tanz 

1+<T 



20 



9^ 



~ 2fl a' dA x 

Ai H xao 1 , - -57-, 

k 1 + cr a£ 



where the factor (w sin i /m) is included for consistency with the convention in our dimensionless variables (which will 
be eventually cancelled). 

MATRICES 

Here we list out the coefficient matrices for the ODEs (see, eqs [55J IH3 US] and 155)) involved in the calculation (in 
the form of A{rj)V{rj) — B(rj)V(?j)) before the reduction procedure in the H4.4I The four ODEs and the five tilde 
variables corresponds to the columns and rows, respectively. We obtain the square matrices by eliminating the fifth 
column with the use of induction equation (see text). Basically the matrices are collections of the background terms in 
the equations and boundary conditions. For the purpose of clarity, we define the following: «t = — v + ft, or = 1 + o 
and wt = oj — (1/L)vt. The "mass matrix", A W) ;: 



where 



b = 



(l+a)'< 



fu T CTT \ 

b u T —xaqS-'/ctt -xao 

ut hxAo 

V u T / 



L , f 2Q 
— a— and h = — 



K \—v/tBia.i 



2f2 tani 



K 1 



a 



(CI) 



(C2) 



The matrix B Wi ; is given by 

/ —ft' — iwx 
2xtoct'/<t t 
(2Cl/K)(il/L)b 

\ 



-<r' i(l/L)(n/2n)a T 
—u' — icuT 1 —xao(1/L) 2 

-(1 + u') -zw T (20/«:)a:Ao/c r T [tani(Z/L) 2 - a'(il/L) 



-(k/2J7) tani 



-iujt 



0\ 




0/ 



(C3) 
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